Attention select de mass (dépendance de mixOmics) prend le pas sur dplyr
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.3 ✓ purrr 0.3.4
## ✓ tibble 3.1.1 ✓ dplyr 1.0.5
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
## Loading required package: lattice
##
## Loaded mixOmics 6.14.1
## Thank you for using mixOmics!
## Tutorials: http://mixomics.org
## Bookdown vignette: https://mixomicsteam.github.io/Bookdown
## Questions, issues: Follow the prompts at http://mixomics.org/contact-us
## Cite us: citation('mixOmics')
##
## Attaching package: 'mixOmics'
## The following object is masked from 'package:purrr':
##
## map
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
##
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
##
## group_rows
chemin <- "../Raw_brut/"
#Echantillons
Sample <- read.csv2(paste0(chemin,"MI_Covariates.csv"))
Corres_var <- read.csv2(paste0(chemin,"MI_Covariates_corresp.csv"))
Corres_ID <- read.csv2(paste0(chemin,"Mb_Info_seq.csv"))
#Métabolites
NMR_CPMG <- read.csv2(paste0(chemin,"NMR_CPMG.csv"))
NMR_NOESY <- read.csv2(paste0(chemin,"NMR_NOESY1D.csv"))
List_cplte_CPMG <- read.csv2(paste0(chemin,"NMR_BUCKET_CPMG_LIST_assigned.csv"))
List_cplte_NOESY <- read.csv2(paste0(chemin,"NMR_BUCKET_NOESY_LIST_assigned.csv"))
List_CPMG <- read.csv2(paste0(chemin,"List_NMR.csv"))
List_NOESY <- read.csv2(paste0(chemin,"List_NOESY.csv"))
## Homogénéisation des noms de variables
##Attention package mass possède un select qui annule celui de dplyr
Sample_ID <- left_join(Sample,Corres_ID,by=c("DonorID"="DonorId")) %>%
relocate("Echantillon", .before = "DonorID") %>%
dplyr::select(-"SeqDepth",-"DonorID")
rownames(Sample_ID) = Sample_ID$Echantillon
#Met la variable AGE en catégorie non hierarchisé et transformation en factor
Sample_ID <- Sample_ID %>%
dplyr::select(-"Echantillon") %>%
mutate(age = factor(AGE, levels = c(1, 2, 3, 4, 5), labels = c("[20-29]", "[30-39]", "[40-49]", "[50-59]", "[60-69]"))) %>%
mutate(SEX = factor(SEX))
Sample_ID$tabac <- factor(Sample_ID$tabac)
CPMG <- left_join(NMR_CPMG,Corres_ID,by=c("DonnorID"="DonorId")) %>%
relocate("Echantillon", .before = "DonnorID") %>%
dplyr::select(-c("DonnorID","Sample_type","VagueManip","Date","Spectre2","Exp","SeqDepth"))
rownames(CPMG) = CPMG$Echantillon
CPMG <- CPMG %>%
dplyr::select(-Echantillon)
NOESY <- left_join(NMR_NOESY,Corres_ID,by=c("DonnorID"="DonorId")) %>%
relocate("Echantillon", .before = "DonnorID") %>%
dplyr::select(-c("DonnorID","Spectre2","Exp","Sample_type","VagueManip","CV_Grps","CV_Grps6","CV_Grps7","CV_Grps8","CV_Grps9","SeqDepth"))
rownames(NOESY) = NOESY$Echantillon
NOESY <- NOESY %>%
dplyr::select(-Echantillon)
# Passage en matrice
Noesy <- as.matrix(NOESY)
Cpmg <- as.matrix(CPMG)
# ?mixOmics
# #selection des métabolites d'interêts. Voir pour transposer et ne selectionner que certaines col avec métabo nommés
# List_CPMG1 <- t(List_CPMG)
# colnames(List_CPMG1) = List_CPMG1[1,]
# List_CPMG2 <- t(List_CPMG1[3,])
# List_CPMG2 <- t(List_CPMG2)tCPMG_vector <- unlist(as.vector(t(CPMG)))
CPMG_vector <- unlist(as.vector(CPMG))
hist(CPMG_vector, breaks = 100)df_row_sum <- apply(CPMG,MARGIN = 1, FUN= mean)
df_col_sum <- apply(CPMG,MARGIN = 2, FUN= mean)
#view(df_sum)
meanr <- mean(df_row_sum)
meanc <- mean(df_col_sum)
minr <- min(df_row_sum)
maxr <- max(df_row_sum)
minc <- min(df_col_sum)
maxc <- max(df_col_sum)
medr <- median(df_row_sum)
medc <- median(df_col_sum)
df <- data.frame(
nb_ech = c(846,202),
moy = c(meanr,meanc),
median = c(medr,medc),
min = c(minr,minc),
max = c(maxr, maxc)
)
rownames(df) <- c("Echantillon","Bucket")
kable(df, align = "c", caption = "Tableau de moyenne des échantillons et buckets") %>%
kable_styling(
#bootstrap_options = "striped",
full_width = F,
position = "center")| nb_ech | moy | median | min | max | |
|---|---|---|---|---|---|
| Echantillon | 846 | 0.0059508 | 0.0059521 | 0.0058323 | 0.0060026 |
| Bucket | 202 | 0.0059508 | 0.0015731 | 0.0000588 | 0.1327695 |
knitr::opts_chunk$set(fig.height = 3.5, fig.width = 8.5)
log10_tCPMG <- log10(t(CPMG) + 1e-8)
log10_CPMG <- log10(CPMG + 1e-8)
# par_ori <- par(no.readonly = TRUE) ## Store the original parameters in order to restore them afterwards
# par(mar = c(4, 6, 5, 1)) ## Adapt the margins of the boxplot in order to have enough space for sample labels
# par(mfrow = c(1,2)) ## Two side-by-side panes
boxplot(t(CPMG),
horizontal = TRUE,
las = 1,
main = "original values",
xlab = "value")Il y a trop de buckets pour que ce soit exploitable. Sachant que seuls 20% sont assignés à des métabolites, il faudrait peut être voir pour ne prendre que ceux là ??
Au vue des boxplot qui sont difficilement exploitables, je suis resté sur la transformation log10(CPMG + 1e-6) Avant transformation, il y avait en sorte 2 groupes : inf et supérieure 0.10 Après transformation cela semble moins marqué mais il reste une “bande” > -1 A mon avis quelque soit la transformation, cela restera. Dois-je continuer pour les centrer ??
Les 2 premiers axes expliquent plus de 62% de la variance exprimée
Quelque soit le groupe étudié (sexe, tabac, âge), sur les axes 1 et 2, les différents points se superposent, il n’est donc pas possible de distinguer un groupe d’un autre. Visualisation montrer en fonction du groupe sexe Il y a un point complétement excentré (18%axe1 et 80% axe2) Homme 60-69ans exfumeur. Doit on garder cet individu (L31R124) ? C’est la plus petite moy (583e-5 et les autres >590e-5) max et min pas si différents des autres
Sur les axes 3 et 4 il est plus facile de distinguer les 2 sexes même s’il y a du chevauchement. Pour les 2 autres il n’est pas posssible de les discrimer. Il y a également 2 points completement excentrés à -20 et -50 axes 3 et 4. On les garde ? Pas discordant dans les autres métadatas
# fviz_pca_ind(res.pca, axes= c(3,4), geom.ind = "point",
# habillage = Sample_ID$age, addEllipses = TRUE)
#
# fviz_pca_ind(res.pca, axes= c(3,4), geom.ind = "point", habillage = Sample_ID$tabac)
La contribution de ces 40 buckets est entre 0.7 et 0.75
fviz_pca_biplot(res.pca, repel = TRUE, axes = c(1,2),
select.var = list(contrib = 40),
geom.ind = "point",
fill.ind = Sample_ID$SEX,
pointshape = 21, pointsize = 2,
palette = "jco",
addaEllipses = TRUE,
col.var = "contrib",
gradient.cols = "BuGn",
legend.title = list(fill = "Sexe", color = "Contrib"))## Warning: ggrepel: 39 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
Pourquoi le graph plotInd de mixOmics semble être symétrique par rapport à celui de factominer L’individu est à-20 sur l’axe1 alors qu’il est à +20 avec ACP de factominer
trans.pca <- pca(log10_CPMG, ncomp=5, center = TRUE, scale = TRUE)
plotIndiv(trans.pca, comp = c(1,2),ind.names = FALSE)Regarde si on peut observer naturellement une clusterisation en fonction d’une catégorie ? Quelque soit la catégorie sexe, age tabac, nous n’observons aucune clusterisation
res.spca.CPMG <- spca(log10_CPMG, scale=TRUE,
ncomp=3, keepX=c(10,10,10))
plotIndiv(res.spca.CPMG, ind.names = FALSE,
group=Sample_ID$SEX,
pch = as.numeric(Sample_ID$SEX),
pch.levels =Sample_ID$SEX,
legend = TRUE)Pourquoi l’axe X variate1 est plus faible que l’axe2. Possibilité de changer manuellement avec comp=c(2,1) X.label
trans.plsda <- plsda(log10_CPMG, Sample_ID$SEX ,ncomp=5, scale = TRUE)
plotIndiv(trans.plsda, ind.names = FALSE, legend=TRUE,
ellipse = TRUE, title="PLS-DA sur CPMG")N’est représenté uniquement les buckets qui contribuent à au moins 70%
plotVar(trans.plsda, var.names = FALSE, cutoff = 0.7,
title = "Correlation circle plot avec cutoff 0.7")Possibilité de faire de la prédiction dans le cas ou les ind seraient en sous groupes, ce qui n’est pas notre cas. ceci est juste un test
trans.testplsda <- plsda(log10_CPMG, Sample_ID$SEX ,ncomp=5, scale = TRUE)
background <- background.predict(trans.testplsda, comp.predicted=2,
dist = "max.dist")
plotIndiv(trans.testplsda, ind.names = FALSE, legend=TRUE,
ellipse = TRUE, title="PLS-DA sur CPMG", background = background)Les variables sélectionnées pour la composante 1 sont exprimées de manière positives pour le sexe1 et négative pour le sexe2
MyResult.plsda2 <- plsda(log10_CPMG, Sample_ID$SEX, ncomp=3)
plotLoadings(MyResult.plsda2, contrib = 'max', method = 'mean')#Ceci est un test pour faire de l'intégration. X = métabo et Y =OTU.
X <- log10_CPMG
Y <- log10(NOESY + 1e-8) ### Il faudrait mettre OTU en Y pour commencer intégration.
MyResult.spls <- spls(X,Y, keepX = c(25, 25), keepY = c(25,25))
plotIndiv(MyResult.spls)## png
## 2
ncomp correspond: The number of components to include in the model for each block (does not necessarily need to take the same value for each block). Je ne sais pas. J’ai mis 3 pour le tabac Ne fonctionne pas pour CPMG
## a Partial Least Squares - Discriminant Analysis is being performed (PLS-DA)
## a Partial Least Squares - Discriminant Analysis is being performed (PLS-DA)
#les 2 ensembles ne fonctionnent pas
#Data_mix_tot <- mixOmics(X=list(Noesy,Cpmg), Y=Sample_ID$tabac, ncomp=3)
#Message :
# Error in Check.entry.wrapper.mint.block(X = X, Y = Y, indY = indY, ncomp = ncomp, :
# Each block of 'X' must have a unique name.
#head(Sample_ID$tabac)
#help(package='mixOmics')